Neutron matter at finite temperature 

L. Tolos^, B. Friman^ and A. Schwenk^ 

^FIAS, J.W. Goethe Universitdt, Ruth-Moufang-Str. 1, 
D-6O438 Frankfurt am Main, Germany 
^GSI, Planckstr. 1, D-64291 Darmstadt, Germany 
^TRIUMF, 4004 Wesbrook Mall, Vancouver, BC, Canada, V6T 2A3 



Abstract 

We calculate the neutron matter equation of state at finite temperature based on 
low-momentum two- and three- nucleon interactions. The free energy is obtained 
from a loop expansion around the Hartree-Fock energy, including contributions from 
normal and anomalous diagrams. We focus on densities below saturation density 
with temperatures T ^ lOMeV and compare our results to the model-independent 
virial equation of state and to variational calculations. Good agreement with the 
virial equation of state is found at low density. We provide simple estimates for the 
theoretical error, important for extrapolations to astrophysical conditions. 



1 Introduction 

The nuclear equation of state plays a central role in astrophysics, for problems 
ranging from the structure of neutron stars neutron star mergers j2] to 
core-collapse supernovae [3)1^ . Astrophysical applications probe the equation 
of state at the extremes of isospin and temperature: The mass of a neutron 
star depends mainly on the equation of state of neutron matter up to densities 
p ~ 4po [5], where po = 0.16 fm~^ is the saturation density of symmetric 
nuclear matter, while supernova explosions are most sensitive to the properties 
of nucleonic matter at subnuclear densities and MeV temperatures [3]. For 
many regimes of interest, the equation of state has to be extrapolated from 
the conditions reached with existing and upcoming experimental facilities. 
Therefore, reliable theoretical input is needed. In this paper, we present a 
study of neutron matter at finite temperature, as part of a program to improve 
the nuclear equation of state input for astrophysics. 
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Conventional nucleon-nucleon (NN) interactions are nonperturbative as a re- 
sult of several sources. First, there is a strong short-range repulsion, which 
leads to bound states of the "flipped" potential AVnn for small, negative A. 
This implies that A = 1 is far outside the radius of convergence. Consequently, 
at least the summation of particle-particle ladder diagrams is required [6] . Sec- 
ond, the tensor force, which is singular at short distances, requires iteration in 
the triplet channels [7f8] . Finally, there are physical bound and nearly-bound 
states in the S-waves, which render the perturbative Born series divergent. Re- 
cently, it was shown that the first two sources of nonperturbative behavior de- 
pend on the choice of NN interaction, and can be removed by evolving nuclear 
forces to low-momentum interactions Viowfc Pll0llllfl2] with cutoffs around 
2fm~^ p^lITi] . An important additional advantage is that the corresponding 
leading-order three-nucleon (3N) interactions from chiral effective field theory 
(EFT) become perturbative in light nuclei for cutoffs A < 2fm~^ [15] . 

At sufficient density (p > O.Olpo [IS] in nuclear matter), Pauli blocking elim- 
inates the shallow bound states, and thus the particle-particle channel be- 
comes perturbative |13] . Consequently, the Hartree-Fock (HF) approximation 
is a good starting point for low-momentum NN and 3N interactions, and per- 
turbation theory (in the sense of a loop expansion) around the HF energy 
becomes tractable. The perturbative character is due to a combination of 
Pauli blocking and an appreciable effective range (see also Ref. [H]). The 3N 
interaction is essential for nuclear matter saturation [13], while the contribu- 
tions to the potential energy remain compatible with EFT power-counting 
estimates. Furthermore, the equation of state becomes significantly less cut- 
off dependent with the inclusion of the dominant second-order contributions. 
In this paper, we extend the investigation of Ref. [13] to neutron matter at 
subsaturation densities, p < po, and generalize the perturbative approach to 
finite temperature. 

Based on the work of Kohn, Luttinger and Ward [TSIITU] , at finite temperature 
the loop expansion around the HF free energy can be realized by the pertur- 
bative expansion of the free energy. In this paper, we include the first-order 
NN and 3N contributions, as well as anomalous and normal second-order dia- 
grams with NN interactions. We defer 3N contributions beyond the HF level 
and higher-order corrections to future work. The pressure, entropy and energy 
are calculated using standard thermodynamic relations. Since low-momentum 
interactions are energy independent, the Matsubara sums can be carried out 
analytically. 

Low- momentum interactions V^owfc and the corresponding 3N forces are defined 
by sharp or smooth regulators with a variable momentum cutoff A. Varying 
the cutoff is a powerful tool to estimate the theoretical errors due to neglected 
higher-order many-body interactions and to assess the completeness of the 
calculations. We use the cutoff dependence to provide simple error estimates. 
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and find that the cutoff dependence is reduced significantly, when second-order 
contributions are included. The possibility of estimating theoretical errors is 
an important step towards reliable extrapolations to astrophysical conditions. 
Finally, we compare our results to the virial equation of state [2011^ and to 
variational calculations [22]. The low-density behavior is in good agreement 
with the virial equation of state. Our results for the energy per particle (see 
Fig. [5]) highlight the importance of a correct finite-temperature treatment of 
second and higher-order correlations. 

This paper is organized as follows. In Sect. [21 we discuss the perturbative 
expansion at finite temperature and give the expressions for the evaluated 
diagrams. Our results for the free energy, pressure, entropy and energy are 
presented in Sect. [31 We conclude and give an outlook in Sect. [H 



2 Loop expansion at finite temperature 

We consider the perturbative expansion of the grand-canonical potential, 

fi(/i,T,V) = -/5 1nZ(/i,T,r), (1) 

where ^(/i, T, V) denotes the partition function of the interacting Fermi sys- 
tem, /i is the chemical potential, f3 = 1/T the inverse temperature and V 
the volume. We include the first-order NN and 3N contributions, fii nn and 
f2i_3N, as well as the second-order anomalous and normal contributions with 
NN interactions, fl2,a. and Q2,n- The grand-canonical potential is then given by 

^ = ^0 + ^1+^2 + ... 

= fio + (^^l,NN + ^^1,3n) + (^^2,a + ^2,n) + • • • , (2) 

where terms of the same order are enclosed in brackets, and is the grand- 
canonical potential of the non-interacting system, 

^ = -2 T / ^ ln(l + e-^(--)) = -2 / ^ . (3) 

V J (27r)3 V J J (27r)3 3m ^ ^ 

Here, = k"^/ (2m) is the free single-particle energy, with m the nucleon mass, 
and Uk = l/[e^*-^'=~^'' + l] is the Fermi-Dirac distribution function. The different 
contributions are depicted diagrammatically in Fig. [H 

The loop expansion around the HF energy is realized by the perturbative 
expansion of the free energy F{N,T,V), which is obtained by a Legendre 
transformation of the grand-canonical potential with respect to the chemical 
potential, 

FiN,T,V) = nif,,T,V)+f,N. (4) 
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Fig. 1. Non-interacting (0); first-order NN (1,NN) and 3N (1,3N); second-order 
anomalous (2, a) and normal (2, n) contributions to the grand-canonical potential 
0,{fi,T,V). The wiggly and double-wiggly lines denote antisymmetrized low-mo- 
mentum NN and 3N interactions, respectively. 



The mean particle number = (A^) is given by 
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In order to invert Eq. ^ for the chemical potential fi{N,T, V), we follow the 
treatment of Kohn and Luttinger [TH] and expand fi to the same order 



(6) 



/i = /iO + /^l + + • • • 



where the particle number is counted as order zero. The lowest order term /iq 
is the chemical potential of a non-interacting system with the same density 
p = N/V as the interacting system. 

Expanding each term on the right-hand side of Eq. around = /xq and 
solving for the chemical potential order by order leads to 
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and /ii 
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and correspondingly for the free energy 



F = Qo{po) + {pi + P2) 
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+ poN + {pi + p2)N + ... . 
Using Eq. ([7j), we find 



F{N)=Fo{N) + ni{p^) + n2{pG 



1 {dnjdpf 



2 d^no/dp^' 



(9) 



Mo 



where Fq{N) = Qq{pq) + pqN is the free energy of the non-interacting sys- 
tem. Since we neglect the contribution of 3N interactions in second-order di- 
agrams, we consistently keep only the first-order NN part fii,NN in the term 
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[dVti/ dij)'^\^Q. Consequently, we have 



F{N) = Fo{N) + fii(/io) + fi2,n(/io) + 



1 (9^VN/M! 

2 92fio/<9/i2 



MO 



• (10) 



According to the Kohn-Luttinger-Ward theorem [TSfTU] , the term in the square 
brackets vanishes at T = for isotropic, normal Fermi systems, since the 
anomalous diagram cancels against the contribution induced in fii,NN by the 
shift of the chemical potential. Thus, the above expansion ensures that the 
T — > limit is correctly reproduced. 

In the T — > limit /iq = ep = A;|/(2m), where kp is the Fermi momentum, and 
the ground state energy of the interacting system is given hj F ^ E = Eq + 
^i{^f) + ^2,n(eF) + . . .. For a momentum-independent contact interaction, the 
square bracket in Eq. ffTOl) vanishes at all temperatures, since in this case the 
HF self-energy is momentum independent. Consequently, the thermodynamic 
potential derived from the free energy Eq. flTU]) corresponds exactly to the 
loop expansion around the HF energy. For finite-range interactions, the HF 
self-energy is momentum dependent, and the cancellation is exact only in the 
zero-temperature limit. At finite temperature, the momentum dependence of 
the HF self-energy is therefore treated perturbatively. 

The pressure, entropy and energy follow from the free energy using standard 
thermodynamic relations. The entropy per particle S/N is given by 



S_ 

N 



d{F/N) 



dT 



d{f/p) 
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where / = F /V is the free-energy density. The chemical potential is given by 
/i = disiF\Ty and the pressure P follows from 



d{F/N) 
'V dN 



d{f/p) 
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dp 
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Finally, the energy per particle is obtained from E/N = F/N + T (S/N). 



2. 1 Hartree-Fock NN and 3N diagrams 



The first-order Viowfc contribution, (1,NN) in Fig. [H is given by 

where the trace is over the spins of the two neutrons and P12 denotes the 
exchange operator for spin and momenta of nucleons 1 and 2. Note that the 
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momentum-conserving delta function is not included in the NN matrix ele- 
ments. 



In neutron matter, the effect of 3N interactions is expected to be smaller 
than in symmetric matter, since the Pauli principle prevents three neutrons 
from interacting in a relative S-state. In the evaluation of the first-order 3N 
diagram, (1,3N) in Fig. [H we follow Ref. [13]. At the HF level only the ci and 
C3 terms of the long-range 27r-exchange part contribute: 



3N 
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138.04 MeV and k 



kj kj. 



As 



where qa = 1-29, /, 
discussed in Ref. [15], we use the Cj constants extracted by the Nijmegen group 
in a partial wave analysis with chiral 27r-exchange [23]: ci = —0.76 GeV~^ and 
C3 = — 4.78GeV~^, where the dominant contribution is due to C3. The low- 
energy constants Cj are within errors consistent with the determination from 
ttN data [21], but at present C3 has a large theoretical uncertainty ~ 25%, 
which is not included in our error bands (see however Fig. [6]). For the 3N 
contribution, we have the regulator [15] 



/r(p, q) = exp 



+ 3gV4^ 
A2 



4n 



(15) 



where p and q are Jacobi momenta. Based on the nuclear matter results of 
Ref. [13] , we expect that the C3 term is repulsive and the dominant part of the 
3N contribution, and that the ci term is small. 



2.2 Second-order anomalous and normal diagrams 



The second-order anomalous contribution, (2, a) in Fig. [H is given by 

V 2T \ l\ "' 7 (27r) 

X (12 I Vio^k (1 - P12) 1 12) (23 I Mowfc (1 - P12) I 23) . (16) 



n Tra,: J 7^ ] Uk, (1 - nfcj nfc3 



We note that in the HF approximation, all tadpole self-energy insertions, 
including the anomalous diagram (2, a), are included in the HF mean- field. 
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Hence, in a loop expansion around the HF solution, the first anomalous 
diagrams are of fourth order and involve two second-order self-energy in- 
sertions in place of the tadpoles in diagram (2, a). As discussed above, the 
present approach is equivalent to a loop expansion around the HF energy for 
a momentum-independent contact interaction. It follows that the contribu- 
tion from the square bracket in Eq. ffTOl) is small, although the anomalous 
diagram (2, a) is significant. The reasons are: First, at zero temperature the 
square bracket vanishes, and therefore the contribution is small at low tem- 
peratures, and second, at finite temperature it is non zero only due to the 
weak momentum dependence of the HF self-energy in neutron matter. 

The second-order normal diagram, (2,n) in Fig. [H reads 



Uk^Uk^ (1 - nk:i){l - UkJ - (1 - nfcj(l - nfca) rik-jUki 



X 



(12|\4owfc(l-A2)|34)'. (17) 



Expanding in partial waves and performing the spin traces, we find 



J2 likSMslVio^kil- Pi2)\^' SM's) 
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TT I 



L J,l,l',S J I J, 

X {k I v,i:i' I k'){k' I v;j:p I k) (i - (-1)^+^+1) (i - (-1)^+^+1) 



X ^J{2l + 1){21' + 1){21 + 1){2T' + 1) (2J + 1)(2J + 1) (-1)'+''+^ 

J lV III L J ) ' 



(/0/'0|L0) (ro/o|Lo 



where 6'kk' is the angle between relative momenta k = (ki — k2)/2 and k' = 
(ka — k4)/2. Keeping only L = in Eq. (fTSl) would result in an angle average 
of the Pauli-blocking operator, but all L ^ 6 are included in our results. 
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Fig. 2. The free energy per particle F/N as a function of density p. The left figure 
gives the first-order NN and 3N contributions with a free single-particle spectrum. 
Second-order anomalous and normal NN contributions are included in the right 
figure. Our results are compared to the virial equation of state (virial) [21] and to 
the variational calculations of Friedman and Pandharipande (FP) [22j . The virial 
curve ends where the fugacity z = e^^'^ = 0.5. 

3 Results 



We compute the different contributions to the free energy using the adaptive 
Monte Carlo integration routine Vegas [25] . Our resulto for the free energy 
per particle are shown in Fig. [2] for temperatures T = 3 MeV, 6 MeV and 
lOMeV, where the low- momentum interaction l^ow/c is obtained from the Ar- 
gonne fig potential [27] for a cutoff A = 2.1 fm~^. The cutoff dependence of 
the free energy can be used to provide lower limits for the theoretical uncer- 
tainty in the calculation, since the result should be cutoff independent when 
all relevant contributions are included. For the T = 6 MeV results, we provide 
error estimates by varying the cutoff over the range A = 1.9 fm~^ (lower curve) 
to A = 2.5 fm~^ (upper curve). The cutoff dependence of the T = 3 MeV and 
10 MeV results here and in the following is of similar size. As expected, the 
error grows with increasing density. Moreover, we observe that the equation 
of state becomes significantly less cutoff dependent with the inclusion of the 
second-order NN contributions. 

In Fig. [2], we also compare our results for the free energy to the model- 
independent virial equation of state [21] and to the variational calculations 

^ We take the opportunity to correct an error in Ref. [26j , where the 3N contribution 
had an incorrect factor in the numerical computation. 
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Fig. 3. The pressure P as a function of density p to first and second order (for details 
see Fig. [2D . Tlie different-temperature FP results are best identified at low density 
by comparison with our results. 

of Friedman and Pandharipande [22] (FP, based on the Argonne W14 and a 
3N potential). The virial expansion provides a benchmark for low densities 
and high temperatures, where the interparticle separation is large compared 
to the thermal wavelength. We find a very good agreement with the virial free 
energy. Including second-order NN contributions to the HF free energy brings 
our results closer to the FP calculations, but this trend is opposite for other 
thermodynamic potentials, see for instance the entropy in Fig. HI 

The FP results are based on zero-temperature Fermi-hypernetted-chain cor- 
relation functions, with the effective mass as a finite-temperature variational 
parameter [28]. We note that the density of states at the Fermi surface is 
underestimated in variational calculations of this type, since the energy de- 
pendence of the self energy is properly accounted for only when correlation 
diagrams are included [29j. In a variational scheme this can be achieved in 
correlated basis perturbation theory [SOfST] . This effect is in part included in 
the induced interaction, which in neutron matter leads to an enhancement of 
the effective mass by ^ 10% (see the RG results for the Fermi liquid parame- 
ter Fi/3 = m* /m — 1 in Fig. 6 of Ref. 22j)- The enhancement of the effective 
mass near the Fermi surface is reflected in an increase of the entropy and the 
specific heat over the variational result at low temperatures [31] (see below). 

The pressure P and the entropy per particle S/N are shown in Figs. [3] and HI 
As for the free energy, we find a very good agreement with the virial equation 
of state at low densities, and the inclusion of second-order contributions signif- 
icantly decreases the cutoff dependence. Our results are similar to the calcu- 
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Fig. 4. The entropy per particle SjN as a function of density p to first and second 
order (for details see Fig. [5]). 
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Fig. 5. The energy per particle EjN as a function of density p to first and second 
order (for details see Fig. [2]) . 

lations of FP for densities p < 0.05 fm~^. For higher densities, we find a larger 
pressure and entropy. We emphasize that the results are based on different 
Hamiltonians, and therefore the comparison has to be taken with care. How- 
ever, the dominant source for the difference in the entropy is likely due to dif- 
ferences in the effective masses, since the entropy density of a low-temperature 
Fermi liquid is proportional to the effective mass, s = m*kpT/3 [33]. As dis- 
cussed above, the variational calculation underestimates the effective mass at 
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0.05 0.1 0.15 



Fig. 6. The energy per particle E/N as a function of density p for T = and 
T = 1.5 MeV. For the T = results, we provide error estimates based on the cutoff 
variation (shaded band) and including the uncertainties in the low-energy constants 
Ci at this level. The upper and lower limits for the range of Cj values are indicated 
by the shaded lines. For comparison we also show the model-independent di-fermion 
EFT results |17j (Schwenk+Pethick) and the results of Akmal et al. [35j . 

the Fermi surface and consequently also the entropy at low temperatures. 

Our results for the energy per particle are presented in Fig. [51 As for the 
free energy, we observe additional binding and a significantly reduced cut- 
off dependence at second order. In contrast to the variational calculation of 
FP [22], the low-density behavior at second order is in good agreement with 
the virial equation of state [21]. This highlights the importance of a correct 
finite-temperature treatment of second and higher-order contributions. Note 
that the error in the virial equation of state (due to the neglected third virial 
coefficient) increases with density. This error is not shown in Fig. O but will be 
discussed in future work on understanding the transition from the perturbative 
to the virial approach. 

A comparison of our low-temperature results to the T = energy per particle 
provides an independent check of our calculations and of the generalized loop 
expansion. In Fig. [6] we show the energy per particle for T = 1.5 MeV for a 
cutoff A = 2.1 fm~^ and the corresponding T = equation of state. The latter 
extends the HF results for neutron matter of Refs. [5^|51] to include 3N forces 
and (normal) second-order NN contributions using an angle-averaged Pauli 
blocking operator (see Ref. [13] for details). Except at low densities, where 
E/N 3/2T, we find that the T = 1.5 MeV energy closely follows the zero 
temperature results. 
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In Fig. [6] we provide an error band for the T = equation of state based 
on the cutoff variation. The width of this band is of the same size as for 
T = 6 MeV in Fig. O At this level of 3N interactions (leading chiral EFT 3N 
in HF), only the 27r-exchange part with low-energy constants q contribute, 
and therefore there are no adjustable 3N parameters. Since the constants 
are cutoff independent, their uncertainties are not fully captured by the cutoff 
variation, and we therefore directly assess how the presently large uncertainties 
in Cj propagate to theoretical uncertainties in the neutron matter equation of 
state at this level. The resulting error estimate in Fig. E] is based on ci = 
-O.QlgJ GeV-^ and Cg = -4.711;^ GeV"^ from Ref. Pj. It is clear that at 
present the theoretical uncertainties in 3N interactions overwhelm the error 
due to an approximate many-body treatment for these densities. 

We can also compare the T = energy per particle at low densities to the 
the model-independent di-fermion EFT results [T7j based directly on the large 
neutron-neutron scattering length and effective range. Our results are consis- 
tent with the di-fermion EFT energy per particle within errors. Finally, the 
results of Akmal et al. [35j (based on the Argonne wis and Urbana IX potential) 
lie within our error band as well (including the q uncertainties). The Urbana 
IX 3N interaction corresponds to the A contribution, = —3.83 GeV~^ [36j, 
and therefore results in less repulsion (with weaker C3). These results show 
that, at present, understanding 3N interactions is a frontier for nuclear mat- 
ter at the extremes. 



4 Conclusions 

This work is part of a program to improve the nuclear equation of state for 
astrophysics. One of the central objectives is to quantify the theoretical un- 
certainties in the microscopic nuclear physics input, and to explore the impact 
on supernovae and neutron stars, for example, through predictions of neutron 
star masses and radii. 

In this first study of neutron matter, we have computed the equation of 
state at subsaturation densities and temperatures T ^ 10 MeV based on low- 
momentum NN and 3N interactions. We have generalized the perturbative 
approach [13] to finite temperature, where the free energy is obtained from a 
loop expansion around the HF energy and the momentum dependence of the 
self-energy is treated perturbatively. Our results include first-order NN and 3N 
contributions, as well as anomalous and normal second-order diagrams with 
NN interactions. The pressure, entropy and energy were then calculated using 
standard thermodynamic relations. While the HF energy is sizable (and non- 
perturbative for finite nuclei), the finite-temperature loop expansion around 
the HF energy seems to be tractable. This is due to a combination of Pauli 
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blocking [13] and an appreciable effective range [T7] . 



The virial expansion provides a model-independent equation of state for nu- 
clear matter at low density and high temperature [201I2T] , and our perturbative 
results meet this benchmark. This is very promising, since it will enable us to 
match the virial equation of state to microscopic calculations based on NN and 
3N interactions at higher densities. The comparison of our results to the virial 
energy per particle highlights the importance of a correct finite-temperature 
treatment of second and higher-order correlations, which are included only 
in an average sense in the variational calculations of Ref. [22] • The correct 
treatment of thermally-excited low-lying states leads to an enhancement of 
the effective mass at the Fermi surface and consequently to an increase in the 
entropy, as shown in Fig. |H 

We have provided simple estimates for the theoretical error by varying the 
cutoff in low-momentum interactions. This is a powerful tool to assess theo- 
retical errors due to neglected higher-order many-body forces and due to an 
approximate many-body treatment. We found that the equation of state be- 
comes significantly less cutoff dependent with the inclusion of second-order 
contributions, and that the cutoff dependence is small for p < 0.1 fm~^. We 
note that the errors of the free energy are correlated between different tem- 
peratures and grow with increasing density. The first observation implies a 
relatively small error in the entropy (obtained by a temperature derivative), 
and consequently similar errors for the energy per particle. The second obser- 
vation explains the relatively large error band for the pressure (obtained from 
a density derivative). Finally, we have shown that the uncertainties due to 
the long-range parts of 3N interactions, the Cj constants, are substantial and 
overwhelm the error bands from the cutoff variation at this level. We conclude 
that understanding 3N forces is a frontier in microscopic calculations of the 
nuclear equation of state, and furthermore that the possibility of estimating 
theoretical uncertainties is an important step towards reliable extrapolations 
to astrophysical conditions. 

Future work will include systematic studies of the range of validity, quantifying 
an expansion parameter for the loop expansion, calculations for asymmetric 
matter, improving the uncertainties at higher densities, how they propagate 
to astrophysical observables, and understanding the transition to the virial 
expansion. 
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